function hs = calcHs_TD(invP,sigmaModel)

A = getA_TD(invP,sigmaModel,invP.disT(1,2));

k = 1./invP.disT(1,2);

hs = cell(invP.nt,1); 
for ii = 1:invP.nt    
   fprintf(1,'Solving time %i of %i....\t',ii,invP.nt);
   if ii == 1;
       prev_hs = invP.h0;
   else
       prev_hs = hs{ii-1};
   end
   rhs = k*invP.mu*prev_hs;
   hs{ii} = A\rhs;
   fprintf(1,'Done\n');
end